
rm(list=ls())
load("TableC3_MSAR_ADL.Rdata")
library(coda)
library(ggplot2)
library(foreign)
class(jags.parsamples) <- "mcmc.list"



epsilon1 <- jags.parsamples[,241][[1]]
epsilon2 <- jags.parsamples[,241][[2]]
epsilon3 <- jags.parsamples[,241][[3]]
combepsilon <- c(epsilon1,epsilon2,epsilon3)
combphi1 <- c(jags.parsamples[,482][[1]],jags.parsamples[,482][[2]],jags.parsamples[,482][[3]])
plot(combepsilon,combphi1)

phi <- rep(NA,7)
for (i in 482:488){
     phi[(i-481)] <- mean(unlist(lapply(jags.parsamples[,i],mean)))
}

delta <- rep(NA,240)
mu <- rep(NA,240)
ypred <- rep(NA,240)
ytilde <- rep(NA,240)

delta <- array(NA,dim=c(240,12000))
mu <- array(NA,dim=c(240,12000))
ypred <- array(NA,dim=c(240,12000))
ytilde <- array(NA,dim=c(240,12000))
for (i in 1:240){
    delta[i,] <- unlist(jags.parsamples[,i])
    mu[i,] <- unlist(jags.parsamples[,(240+1+i)])
    ypred[i,] <- unlist(jags.parsamples[,(490+i)])
    ytilde[i,] <- (mu[i,] + unlist(jags.parsamples[,483])*mean(statadata$party) + unlist(jags.parsamples[,484])*mean(statadata$partyics) + unlist(jags.parsamples[,485])*mean(statadata$partyapprove) + unlist(jags.parsamples[,486])*mean(statadata$lagparty) + unlist(jags.parsamples[,487])*mean(statadata$lagpartyics) + unlist(jags.parsamples[,488])*mean(statadata$lagpartyapprove))/(1-unlist(jags.parsamples[,482]))
}

meandelta <- apply(delta,c(1),mean)
meanytilde <- apply(ytilde,c(1),mean)
meanmu <- apply(mu,c(1),mean)
meanypred <- apply(ypred,c(1),mean)

meanytilde2 <- (meanmu + phi[2]*mean(statadata$party) + phi[3]*mean(statadata$partyics) + phi[4]*mean(statadata$partyapprove) + phi[5]*mean(statadata$lagparty) + phi[6]*mean(statadata$lagpartyics) + phi[7]*mean(statadata$lagpartyapprove))/(1-phi[1])

mean(meanytilde)
mean(meanytilde2)
mean(statadata$adjusted_gallup)

## R-squared measures
sum((meanytilde-mean(statadata$adjusted_gallup))^2)/sum((statadata$adjusted_gallup-mean(statadata$adjusted_gallup))^2)
sum((meanytilde2-mean(statadata$adjusted_gallup))^2)/sum((statadata$adjusted_gallup-mean(statadata$adjusted_gallup))^2)
sum((meanypred-mean(statadata$adjusted_gallup))^2)/sum((statadata$adjusted_gallup-mean(statadata$adjusted_gallup))^2)
res <- statadata$adjusted_gallup - ypred
sqrt(mean(res^2))




statadata$mu <- meanmu
statadata$ytilde <- meanytilde2
statadata$delta <- meandelta
statadata$ypred <- meanypred


mudiff <- mu[2:240,]-mu[1:239,]
mudiff <- mudiff*100
sqmudiff <- mudiff^2
absmudiff <- abs(mudiff)

futurevar <- array(NA,dim=c(240,12000))
futureabsdev <- array(NA,dim=c(240,12000))

for (i in 1:238){
    if (i<233){
        futureabsdev[i+1,] <- apply(absmudiff[i:(i+7),],c(2),sum)
        futurevar[i+1,] <- apply(sqmudiff[i:(i+7),],c(2),sum)/8
    }    else {
        futureabsdev[i+1,] <- apply(absmudiff[i:239,],c(2),sum)
        futurevar[i+1,] <- apply(sqmudiff[i:239,],c(2),sum)/(240-i)
    }
}

futurevar[1,] <- apply(sqmudiff[1:6,],c(2),sum)/7
futureabsdev[1,] <- apply(absmudiff[1:6,],c(2),sum)
futurevar[240,] <- sqmudiff[239,]
futureabsdev[240,] <- absmudiff[239,]

avgfvar <- apply(futurevar,c(1),mean)
avgabsdev <- apply(futureabsdev,c(1),mean)

avgfvartau <-  1/apply(futurevar,c(1),var,na.rm=TRUE)
lbfvar <- apply(futurevar,c(1),quantile,probs=c(.05),na.rm=TRUE)
ubfvar <- apply(futurevar,c(1),quantile,probs=c(.95),na.rm=TRUE)

futuresd <- sqrt(futurevar)
avgfsd <-  apply(futuresd,c(1),mean)
avgfsdtau <-  1/apply(futuresd,c(1),var,na.rm=TRUE)
lb5fsd <- apply(futuresd,c(1),quantile,probs=c(.05),na.rm=TRUE)
lb10fsd <- apply(futuresd,c(1),quantile,probs=c(.10),na.rm=TRUE)
ubfsd <- apply(futuresd,c(1),quantile,probs=c(.95),na.rm=TRUE)

plot(avgfsd,type="l",ylim=c(-0.5,2))
lines(lb5fsd,lty=3)
lines(lb10fsd,lty=3)

statadata$avgfvar <- avgfvar
statadata$avgfvartau <- avgfvartau
statadata$lbfvar <- lbfvar
statadata$avgfsd <- avgfsd
statadata$avgfsdtau <- avgfsdtau
statadata$lb5fsd <- lb5fsd
statadata$lb10fsd <- lb10fsd


combdelta <- rbind(delta,ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0))

futureshift <- array(NA,dim=c(240,12000))

for (i in 1:240){
    futureshift[i,] <- ifelse(apply(combdelta[i:(i+7),],c(2),sum)>0,1,0)
}

avgfutureshift <- apply(futureshift,c(1),mean)
statadata$avgfutureshift <- avgfutureshift


write.dta(statadata,"TableD6_MSAR_ADL_updatingestimates.dta")

rm(list=ls())
                   

